#!/usr/bin/env python3
"""
糖尿病病情严重程度预测 - Auto-sklearn回归建模
==================================================

这是一个完整的医疗AI回归预测示例，展示如何使用Auto-sklearn
预测糖尿病患者的疾病严重程度量化值。

功能特色：
- 全自动机器学习回归建模
- 医疗专用评估指标
- 临床解释性分析
- 模型集成和优化
"""

import os
import sys
import warnings
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import datasets, model_selection, metrics
from sklearn.preprocessing import StandardScaler
import autosklearn.regression
from datetime import datetime
import joblib

# 设置警告过滤和显示
warnings.filterwarnings('ignore')
plt.style.use('seaborn-v0_8')

class DiabetesRegressionPredictor:
    """糖尿病严重程度回归预测器"""
    
    def __init__(self, time_limit=180, per_run_limit=30):
        """
        初始化预测器
        
        Args:
            time_limit: Auto-sklearn总运行时间（秒）
            per_run_limit: 单个模型最大运行时间（秒）
        """
        self.time_limit = time_limit
        self.per_run_limit = per_run_limit
        self.automl = None
        self.scaler = None
        self.feature_names = None
        self.results = {}
        
    def load_and_prepare_data(self):
        """加载和预处理糖尿病数据集"""
        print("正在加载糖尿病数据集...")
        
        # 加载数据
        diabetes = datasets.load_diabetes()
        X, y = diabetes.data, diabetes.target
        
        # 特征名称
        self.feature_names = [
            'age', 'sex', 'bmi', 'bp', 's1', 's2', 's3', 's4', 's5', 's6'
        ]
        
        # 创建DataFrame以便分析
        df = pd.DataFrame(X, columns=self.feature_names)
        df['target'] = y
        
        print(f"数据集形状: {X.shape}")
        print(f"目标变量范围: {y.min():.1f} - {y.max():.1f}")
        print(f"目标变量均值: {y.mean():.1f} ± {y.std():.1f}")
        
        # 训练集和测试集分割
        X_train, X_test, y_train, y_test = model_selection.train_test_split(
            X, y, test_size=0.25, random_state=42, stratify=None
        )
        
        print(f"训练集大小: {X_train.shape[0]} 样本")
        print(f"测试集大小: {X_test.shape[0]} 样本")
        
        return X_train, X_test, y_train, y_test, df
    
    def create_medical_regressor(self):
        """创建医疗专用Auto-sklearn回归器"""
        print("创建Auto-sklearn回归器...")
        
        # 配置AutoML回归器
        self.automl = autosklearn.regression.AutoSklearnRegressor(
            time_left_for_this_task=self.time_limit,
            per_run_time_limit=self.per_run_limit,
            initial_configurations_via_metalearning=25,
            smac_scenario_args={'runcount_limit': 50},
            resampling_strategy='cv',
            resampling_strategy_arguments={'folds': 5},
            tmp_folder="/tmp/autosklearn_diabetes_regression",
            seed=42,
            memory_limit=3072,
            include={
                'regressor': [
                    'random_forest', 'extra_trees', 'gradient_boosting',
                    'ridge_regression', 'lasso_regression', 'elastic_net',
                    'support_vector_regression', 'k_nearest_neighbors'
                ]
            }
        )
        
        return self.automl
    
    def train_model(self, X_train, y_train):
        """训练Auto-sklearn模型"""
        print(f"开始训练模型 (时间限制: {self.time_limit}秒)...")
        print("=" * 50)
        
        start_time = datetime.now()
        
        # 训练模型
        self.automl.fit(
            X_train, y_train,
            dataset_name="diabetes_regression"
        )
        
        end_time = datetime.now()
        training_time = (end_time - start_time).total_seconds()
        
        print(f"模型训练完成，耗时: {training_time:.1f}秒")
        print("=" * 50)
        
        return training_time
    
    def calculate_medical_metrics(self, y_true, y_pred):
        """计算医疗回归专用评估指标"""
        metrics_dict = {}
        
        # 基础回归指标
        metrics_dict['r2_score'] = metrics.r2_score(y_true, y_pred)
        metrics_dict['mse'] = metrics.mean_squared_error(y_true, y_pred)
        metrics_dict['rmse'] = np.sqrt(metrics_dict['mse'])
        metrics_dict['mae'] = metrics.mean_absolute_error(y_true, y_pred)
        metrics_dict['mape'] = np.mean(np.abs((y_true - y_pred) / y_true)) * 100
        
        # 医疗相关指标
        # 相关系数
        correlation = np.corrcoef(y_true, y_pred)[0, 1]
        metrics_dict['correlation'] = correlation
        
        # 解释方差得分
        metrics_dict['explained_variance'] = metrics.explained_variance_score(y_true, y_pred)
        
        # 临床可接受性评估 (误差在20%以内的比例)
        relative_errors = np.abs((y_true - y_pred) / y_true)
        metrics_dict['clinical_acceptable_rate'] = np.mean(relative_errors <= 0.2) * 100
        
        # 预测稳定性 (标准差)
        residuals = y_true - y_pred
        metrics_dict['prediction_stability'] = np.std(residuals)
        
        return metrics_dict
    
    def evaluate_model(self, X_train, X_test, y_train, y_test):
        """全面评估模型性能"""
        print("正在评估模型性能...")
        
        # 生成预测
        train_pred = self.automl.predict(X_train)
        test_pred = self.automl.predict(X_test)
        
        # 计算评估指标
        train_metrics = self.calculate_medical_metrics(y_train, train_pred)
        test_metrics = self.calculate_medical_metrics(y_test, test_pred)
        
        # 保存结果
        self.results.update({
            'train_predictions': train_pred,
            'test_predictions': test_pred,
            'train_metrics': train_metrics,
            'test_metrics': test_metrics
        })
        
        # 显示评估结果
        print("\n📊 模型性能评估结果")
        print("=" * 60)
        print(f"{'指标':<25} {'训练集':<15} {'测试集':<15}")
        print("-" * 60)
        print(f"{'R² 决定系数':<25} {train_metrics['r2_score']:<15.4f} {test_metrics['r2_score']:<15.4f}")
        print(f"{'RMSE':<25} {train_metrics['rmse']:<15.2f} {test_metrics['rmse']:<15.2f}")
        print(f"{'MAE':<25} {train_metrics['mae']:<15.2f} {test_metrics['mae']:<15.2f}")
        print(f"{'MAPE (%)':<25} {train_metrics['mape']:<15.1f} {test_metrics['mape']:<15.1f}")
        print(f"{'相关系数':<25} {train_metrics['correlation']:<15.4f} {test_metrics['correlation']:<15.4f}")
        print(f"{'临床可接受率 (%)':<25} {train_metrics['clinical_acceptable_rate']:<15.1f} {test_metrics['clinical_acceptable_rate']:<15.1f}")
        
        # 模型性能诊断
        self._diagnose_performance(test_metrics)
        
        return train_metrics, test_metrics
    
    def _diagnose_performance(self, metrics):
        """诊断模型性能并提供建议"""
        print("\n🔍 模型性能诊断")
        print("=" * 40)
        
        r2 = metrics['r2_score']
        rmse = metrics['rmse']
        clinical_rate = metrics['clinical_acceptable_rate']
        
        if r2 >= 0.7:
            print("✅ 优秀：模型解释了70%以上的方差")
        elif r2 >= 0.5:
            print("✅ 良好：模型解释了50%以上的方差")
        elif r2 >= 0.3:
            print("⚠️  一般：模型性能需要优化")
        else:
            print("❌ 差：模型性能不佳，需要重新设计")
        
        if clinical_rate >= 80:
            print("✅ 临床可接受：80%以上预测在可接受范围内")
        elif clinical_rate >= 60:
            print("⚠️  临床谨慎：仅60%以上预测可接受")
        else:
            print("❌ 临床不可接受：预测准确性不足")
        
        if rmse <= 30:
            print("✅ 预测精度：RMSE在合理范围内")
        else:
            print("⚠️  预测精度：RMSE偏高，需要改进")
    
    def analyze_model_insights(self):
        """分析模型洞察和可解释性"""
        print("\n🔬 模型分析和洞察")
        print("=" * 50)
        
        # 显示最佳模型
        print("🏆 最优模型排行榜:")
        print(self.automl.leaderboard())
        
        print("\n🔧 最终集成模型构成:")
        print(self.automl.show_models())
        
        # 模型统计信息
        print(f"\n📈 模型搜索统计:")
        print(f"- 评估的配置数量: {len(self.automl.cv_results_)}")
        print(f"- 最佳交叉验证分数: {self.automl.cv_results_['mean_test_score'].max():.4f}")
    
    def create_visualizations(self, X_test, y_test):
        """创建可视化图表"""
        print("正在生成可视化图表...")
        
        test_pred = self.results['test_predictions']
        
        # 创建图表
        fig, axes = plt.subplots(2, 2, figsize=(15, 12))
        fig.suptitle('糖尿病严重程度预测分析', fontsize=16, fontweight='bold')
        
        # 1. 预测 vs 真实值散点图
        ax1 = axes[0, 0]
        ax1.scatter(y_test, test_pred, alpha=0.6, color='steelblue')
        ax1.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], 'r--', lw=2)
        ax1.set_xlabel('真实值')
        ax1.set_ylabel('预测值')
        ax1.set_title('预测值 vs 真实值')
        ax1.grid(True, alpha=0.3)
        
        # 添加R²注释
        r2 = self.results['test_metrics']['r2_score']
        ax1.text(0.05, 0.95, f'R² = {r2:.3f}', transform=ax1.transAxes, 
                bbox=dict(boxstyle="round", facecolor='wheat'), fontsize=12)
        
        # 2. 残差图
        ax2 = axes[0, 1]
        residuals = y_test - test_pred
        ax2.scatter(test_pred, residuals, alpha=0.6, color='green')
        ax2.axhline(y=0, color='r', linestyle='--')
        ax2.set_xlabel('预测值')
        ax2.set_ylabel('残差')
        ax2.set_title('残差分析')
        ax2.grid(True, alpha=0.3)
        
        # 3. 残差分布直方图
        ax3 = axes[1, 0]
        ax3.hist(residuals, bins=20, alpha=0.7, color='orange', edgecolor='black')
        ax3.set_xlabel('残差')
        ax3.set_ylabel('频次')
        ax3.set_title('残差分布')
        ax3.axvline(x=0, color='r', linestyle='--')
        ax3.grid(True, alpha=0.3)
        
        # 4. 误差百分比分布
        ax4 = axes[1, 1]
        relative_errors = np.abs((y_test - test_pred) / y_test) * 100
        ax4.hist(relative_errors, bins=20, alpha=0.7, color='purple', edgecolor='black')
        ax4.set_xlabel('相对误差 (%)')
        ax4.set_ylabel('频次')
        ax4.set_title('预测误差分布')
        ax4.axvline(x=20, color='r', linestyle='--', label='临床可接受阈值')
        ax4.legend()
        ax4.grid(True, alpha=0.3)
        
        plt.tight_layout()
        
        # 保存图表
        os.makedirs('../results', exist_ok=True)
        plt.savefig('../results/regression_analysis.png', dpi=300, bbox_inches='tight')
        print("✅ 图表已保存至 ../results/regression_analysis.png")
        
        return fig
    
    def save_model_and_results(self):
        """保存模型和结果"""
        print("正在保存模型和结果...")
        
        os.makedirs('../results', exist_ok=True)
        
        # 保存Auto-sklearn模型
        model_path = '../results/diabetes_automl_model.pkl'
        joblib.dump(self.automl, model_path)
        
        # 保存评估结果
        results_path = '../results/evaluation_results.pkl'
        joblib.dump(self.results, results_path)
        
        print(f"✅ 模型已保存至: {model_path}")
        print(f"✅ 结果已保存至: {results_path}")
    
    def generate_clinical_report(self):
        """生成临床应用报告"""
        print("\n📋 临床应用报告")
        print("=" * 50)
        
        test_metrics = self.results['test_metrics']
        
        print("🎯 预测性能摘要:")
        print(f"- 模型解释方差: {test_metrics['r2_score']:.1%}")
        print(f"- 平均预测误差: ±{test_metrics['mae']:.1f} 单位")
        print(f"- 临床可接受率: {test_metrics['clinical_acceptable_rate']:.1f}%")
        
        print("\n💡 临床应用建议:")
        if test_metrics['r2_score'] >= 0.5:
            print("✅ 推荐用于临床辅助诊断")
            print("✅ 可用于患者风险分层")
            print("✅ 建议结合医生专业判断使用")
        else:
            print("⚠️  建议进一步优化后再用于临床")
            print("⚠️  需要更多数据或特征工程")
        
        print("\n🔄 持续改进方向:")
        print("- 收集更多临床特征数据")
        print("- 集成时间序列信息")
        print("- 考虑患者个体差异")
        print("- 建立模型更新机制")

def main():
    """主函数"""
    print("🏥 糖尿病严重程度预测系统")
    print("=" * 50)
    
    # 创建预测器
    predictor = DiabetesRegressionPredictor(time_limit=180, per_run_limit=30)
    
    try:
        # 1. 数据加载和预处理
        X_train, X_test, y_train, y_test, df = predictor.load_and_prepare_data()
        
        # 2. 创建和训练模型
        predictor.create_medical_regressor()
        training_time = predictor.train_model(X_train, y_train)
        
        # 3. 模型评估
        train_metrics, test_metrics = predictor.evaluate_model(X_train, X_test, y_train, y_test)
        
        # 4. 模型分析
        predictor.analyze_model_insights()
        
        # 5. 可视化
        predictor.create_visualizations(X_test, y_test)
        
        # 6. 保存结果
        predictor.save_model_and_results()
        
        # 7. 生成临床报告
        predictor.generate_clinical_report()
        
        print(f"\n🎉 分析完成！总耗时: {training_time:.1f}秒")
        print("📂 结果文件保存在 ../results/ 目录中")
        
    except Exception as e:
        print(f"❌ 运行出错: {str(e)}")
        raise

if __name__ == "__main__":
    main()